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Abstract. Hydrodynamic unstratified keplerian flows are known to be linearly stable at all Reynolds numbers, but may never- 
theless become turbulent through nonlinear mechanisms. However, in the last ten years, conflicting points of view have appeared 
on this issue. We have revisited the problem through numerical simulations in the shearing sheet limit. It turns out that the effect 
of the Coriolis force in stabilizing the flow depends on whether the flow is cyclonic (cooperating shear and rotation vorticities) 
or anticyclonic (competing shear and rotation vorticities); keplerian flows are anticyclonic. We have obtained the following 
results: 

i/ The Coriolis force does not quench turbulence in subcritical flows; however, turbulence is more efficient, and much more 
easily found, in cyclonic flows than in anticyclonic ones. 

ii/ The Reynolds number/rotation/resolution relation has been quantified in this problem. In particular we find that the resolution 
demand, when moving away from the marginal stability boundary, is much more severe for anticyclonic flows than for cyclonic 
ones. Presently available computer resources do not allow numerical codes to reach the keplerian regime, 
hi/ The efficiency of turbulent transport is directly correlated to the Reynolds number of transition to turbulence Rg, in such 
a way that the Shakura-Sunyaev parameter a ~ 1/Rg. This correlation is nearly independent of the flow cyclonicity. The 
correlation is expected on the basis of generic physical arguments. 

iv/ Even the most optimistic extrapolations of our numerical data show that subcritical turbulent transport would be too inef- 
ficient in keplerian flows by several orders of magnitude for astrophysical purposes. Vertical boundary conditions may play a 
role in this issue although no significant effect was found in our preliminary tests. 

v/ Our results suggest that the data obtained for keplerian-like flows in a Taylor-Couette settings are largely affected by sec- 
ondary flows, such as Ekman circulation. 
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H | 1. Introduction Hawley and their co llaborators jBalbus & Hawlevl 1 199 it 

■ - - > lHawlev et alJI 19951 see lBalbusll2003l for a recent review - ). The 

The question of the existence and physical origin of turbu- . , , . . . • , . . • , ■ 

M turbulent transport induced by this instability is by now charac- 

lence in accretion disks has been lively debated for a number of . . , . juu n j 

3 tenzed in a number of instances, and has been called upon even 
decades. Generally speaking, there are a priori two basic ways , , c .. c ,. , . . . , . ., 

3 r b r 3 when only some fraction of the disk is ionized, as in the mid- 

in which an accretion disk can become turbulent. In the first , e vc< -. • .• , tU . . Tp, H 

plane region of YSOs inner disks — the dead-zone ( Gammie 

way, some linear instability is present in the flow, and its non- m& plemin2 & Ston£ 200 | Howeyeri the reduced efficiency 

linear development eventually drives turbulence. In the second , , . ... ,, U1 • . 

F 3 of the transport in this case, as well as the possible existence 

one, the flow is linearly stable, and undergoes a direct laminar- c ,. , , . , , , uim , , , 

3 . of disks which may not support MHD phenomena at all, has 

turbulent transition once a certain threshold in Reynolds num- . , c . . . ■ iujj 

3 prompted some upsurge of interest in purely hydrodynamic 

ber is reached. The first type of transition to turbulence is called • . A , , , .. . ri . .... , , , 

JF instabilities. A local, barochmc-hke instability has been ob- 

supercritical, and the second, (globally) subcritical. , . , • w - , I^TT s> t> a u ■ I Jonml 

F ' ' V6 I — n — 1 1 served in global simulations by Klahr & Bodenheimei (2003). 

Global instabilities (such as the| Papaloizou & P ringle 198J Local stabilitv analyzes (KlahJ l2004 Llohnson & Gammie 

instability) seem unpromising to drive turbulence <BlaeJl98j EoQ53t find transient instability in this context, but shear- 

| Hawley || l99l|) . As for local instabilities, an astrophysically ing box simu l a tions indicate that this does not drive tur- 

important example of supercritical transition is provided by bulence J Johnson & G ammiel Ebo5b). Urmn (2003) discusses 

the magneto-rotational instability (MRI) which has been ex- an instability related t0 vertic al shear and heat transpor t 

tensively studied followi ng the pioneering work of Balbus, of ±e Goldrei ch-Schubert type fctoldreich & Sch^fl%7i) : 

Send offprint requests to: G. Lesur however, this instability produces only a rather weak radial 
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transpo rt dArlt & Urpinll2004. More recen tly, iDubrulle etaD 
J2005bl) and lShalvbkov & Ruedigerl §005) have discussed an 
instability arising when both the fluid differential rotation and 
vertical stratification are stabilizing according to the H0iland 
criterion. However, it seems that this instability is connected 
to the presence of walls, and is dynamically important only 
when the inter-wall distance is small enough for a resonant- 
like interaction to take place 1 (Satomura 1981), otherwise dis- 
turbances are confined to the near boundary zone; a related 
result has recen tly been found in the astrophysics literature 
JUmurhanl2005l) . Earlier analytic and numerical investigations 
have shown this instability to be absent in local disk mod- 
els JGoodman & B albus 20oH [Brandenburg & Dintransll200ll 



This conclusion was in turn questioned by Richa rd & Zahnl 



Rudi ger et alJl2002h . Note finally that vertical convection in a 
stratified disk can in principle also drive turbulence; however, it 
indu ces inwards transport ins tead of the required outwards one 
JCaboll996HStone & Balbusll99fj|) . Therefore, no local insta- 
bility has yet been found in the hydrodynamic regime, which 
would explain the turbulent transport taking place in accretion 
disks. 

Subcritical transition to turbulence is the subject of the 
present work. The non-rotating plane Couette flow provides 
a classical (and to date the best understood) example of a 
system undergoing a subcritical transition. Although the na- 
ture and mechanism of the transition remained elusive for 
decades, it ha s been identified in the recent years, in laboratory 
experiments ( Daviaud et al. 1992; Dauchot & Daviaud 1995a; 
iDauchot & Daviaud|l995 bHBottin et alJl997|)!ri umerica l sim- 
ulatio ns lHa^rultonetal*|T 995t see also l Schmiegel & Eckh ardt 
1 19971 and lEckhardt & Me" smann 1999), and theoretical ana- 
lyzes (in particular Waleffe 1997; Waleffe 2003). Earlier in- 
vestigations of the problem have focused on the role of nonlin- 
ear instabilities in subcritical shear flo ws, based on Landau- 
like toy-models on the one hand (e.g.. iDrazin & Reidlfl98ll 
and references therein), and analysis of the linear stability 
of finite a mplitude defect s in the flow profile on the o ther 
llLerner & Knoblochlll988i IDubrulle & Zahnlll99l[ IPubrullel 
Il993l) : unfortunately, such analyzes yield little information on 
the existence and location of the turbulent state in parameter 
space and on the turbulent transport efficiency, unless further 
ad hoc assumptions are made. 

In any case, on the basis of the empirically observed sub- 
critical transition in laboratory flows, it was suggested that a 
similar process is relevant in accretion disks (Sh akura et alJ 
Il978h . in spite of their very different prevailing physical condi- 
tions. This suggestion was tested and challenged in a seri es of 
numerical sim ulations performed by Bal bus et alJ (|1996) and 
Hawl ev et al 1 dl999i) . in the shearing sheet limit. Transition to 
turbulence was not found in these simulations for keplerian- 
like flows. The simulations were performed with two different 
finite difference codes (a PPM type code, and the ZEUS code), 
up to a resolution of 256 3 . These two works concluded that a 
stabilizing Coriolis force prevents the existence of turbulence 
in the simulated flows, except in the immediate vicinity of the 
linear marginal stability limits. 



( 1999J), o n the basis of the Taylor-Couette experiments per- 
formed bv lWendtH 1933b and lTavloJ (fT936). These experimen- 
tal results display a subcritical transition to turbulence in pres- 
ence of a stabilizing Coriolis force. Also, new sets of experi- 
ments have been carried out in order to bring the experimen- 
tal conditions closer to the ones prevailing in a keplerian flow. 
Namely, a Taylor-Couette apparatus was used in conditions 
of radially decreasing angular velocity and radially increasing 
specific angular momentum. Turbulence was again found for 
high enough Reynolds numbers dRicharcfeOOlt iRichard et alJ 
2001) but the results are not unambiguous, as the potential role 
of secondary flows induced by the boundary conditions in the 
experiments, such as Ekmann's circulation, is unclear, in spite 
of the attention devoted to this point in the experiments. In any 
case, a subcritical transition is also found in all experiments 
of shear flows on which a linearly stabilizing Coriolis force is 
su perimposed (lLongaretti & Pauchoj2005b . 

iLongarettH J2002I) has argued from a phenomenological 
analysis that the lack of turbulence in the simulations per- 
formed to date was due to a lack of resolution, as the Coriolis 
force may increase the range of scales that need to be resolved 
for a subcritical turbulent transition to show up. On the other 
hand, on the basis of a newly deve loped Reynolds stress closure 
scheme (IO gilviel2003b . lCjaraud & Qgilviel ll2005b find that ke- 
plerian flows may or may not be turbulent depending on the pa- 
rameters of the scheme. For their favored choice of parameters, 
unbounded keplerian flows are not turbulent, on the contrary to 
linearly stable, wall-bounded Taylor-Couette flows. 

The recent astrophysical literature on the problem of sub- 
critical transition has also focuse d on the concept of tran 
sient growth in keplerian flows ( Chagelishvili et al] 20Tj2 
)03ll^ckoll2004IUmurhan&Regevl 



1 We thank Stephane Le Dizes for bringing this point to our atten- 
tion. 



Tevzadze et al. 2003; Yecko 2004; Umurhan & Regev 2004; 
Mukhopadhv av et al.ll2004l (Afshordi et al. 2004} . Pue to the 
nonnormal character of the Navier-Stokes equation, linear 
modes can transiently be strongly amplified in shear flows, al- 
though on the long run they must viscously decay. It has been 
argued that this transient growth can be relevant to astrophys- 
ical disks in two different ways. First, 3P turbulence (or an 
external forcing) can couple to large scale 2P structures; the 
(statistical) amplitude of these structures can be large, under 
the combined action of this coupling, of transient growth and 
of viscous decay, and these 2P structures may contribute t o 
the overall transport in the disk Jloannou & Kakourisll200ll) . 
Secondly, a large transient growth has been invoked in the by- 
pass scenario of transition to turbulence, which involves an 
interplay between nonnormality and nonlin e arity (se e , e.g. , 
ICirossmaiJEoOoT: iBrosa & Grossmannlll999l) . IWaleffd J 19951) 
has emphasized the key role played by nonlinear interactions in 
the context of the recently identified turbulent self-sustaining 
process of non-rota ting plane Couette flows ^Hamilton et alJ 
1995; Waleffe 1997). Even though transient growth explains 
the strong modulations of the streamwise velocity from rel- 
atively weak streamwise rolls involved in this self-sustaining 
mechanism, the existence and properties of the turbulent basin 
of attraction for the full nonlinear dynamics are apparently 
poorly constrained by the nonnormal linear problem charac- 
teristics. 
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Our present understanding of the possible existence of a 
dynamically significant subcritical turbulent transition in ac- 
cretion disks is unsatisfying in several respects, calling for a 
reinvestigation of the problem. On the one hand, the relevance 
of the available laboratory experiments to accretion disk tur- 
bulence is at best unclear, as will be shown in the course of 
the p resent work (for a different opinion, see iHersant et alJ 
120051) . On the other hand, the absence of subcritical turbulence 
in the shearing she et loc al model of accretio n disks used by 
Balbus et al. ( 1996J and lHawlev et all ll 19991) may be an ef- 
fect of various numerical limitations, namely, algorithm choice, 
limited resolution, nature of the boundary conditions, imposed 
aspect ratio and initial conditions of the simulations. Of these 
options, only the first two have been partially addressed in 
these previous investigations, leading to questions concerning 
the "effective Reynolds number" of the performed simulations 
— an ill-defined process-dependent concept, that we shall clar- 
ify in t he context of the p resent problem. Following the sugges- 
tion of lLongarettil d2002h . the primary aim of the present work 
is to investigate in a more systematic way, through numerical 
simulations of plane parallel, rotating shear flows, the effects 
of finite resolution on the results. The effects of the other fac- 
tors listed above are also somewhat explored, but to a lesser 
extent. Both cyclonic and anticyclonic rotation are considered; 
although cyclonic rotation is not relevant to accretion disks, it 
turns out that cyclonic flows behave very differently from an- 
ticyclonic ones, opening some interesting perspective into the 
nature of the problem. 

This paper is organized as follows. Section l2~Tl collects the 
background material relevant to the problem. First, the form of 
the equations solved is provided, and the global energy budget 
recalled, before discussing linear stability limits. The section is 
concluded by a summary of the effect of a stabilizing rotation in 
shear flows as characterized by the available laboratory experi- 
ments. The next section presents the various codes used in this 
work, and the numerical results obtained with them. Section^] 
discusses various aspects of our numerical results, most notably 
the role of resolution and boundary conditions on the numerical 
side, the role of the Coriolis force, the underlying phenomeno- 
logical picture, and the astrophysical implications, on the phys- 
ical side. A summary is provided in section [5] along with an 
outlook on the question of turbulence in accretion disks. 

2. Rotating plane shear flows: a summary 

The present investigation is concerned with the nonlinear in- 
stability of laminar flows characterized by a uniform shear, in 
the presence of a uniform global rotation. The direction of the 
flow is identified with the x axis (streamwise direction), and 
the direction of the shear with the y axis (shearwise direction); 
rotation is applied along the z axis (spanwise direction). The 
laminar flow u^, is invariant in the streamwise and spanwise di- 
rections (in particular, the vertical stratification expected in a 
real disk is ignored): u^ = U(y)e x . 

Such a flow can be used to numerically model either a local 
portion of an accretion disk, or experiments on rotating plane 
Couette flows, depending on the nature of the applied boundary 
condition in the shearwise direction (in practice, either rigid or 



shearing sheet; see next section). The configuration is repre- 
sented on Fig.^ 
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Fig. 1. Sketch of the configuration of rotating plane shear flows. 



2.1. Equations of motion 

The most useful form of the Navier-Stokes equation, for our 
present purpose, is obtained by separating the laminar flow Ul 
and the deviation from laminar w in the total velocity u in the 
rotating frame, leading to 



<9w <9w 

— + w • Vw = S ■ y— + (2Q. + S)w y e x - 2Qw t e y 
ot ox 



VSk 



vAw, 



(1) 



where the gradient terms balancing the laminar flow Coriolis 
force has been subtracted out to form the effective generalized 
pressure 5n (which therefore absorbs the equilibrium centrifu- 
gal, gravitational and/or pressure force term, depending on the 
considered equilibrium problem); Q is the flow rotation veloc- 
ity in an inertial frame, and S = -dU/dy is the shear. The 
convention adopted here is that the sign of S is chosen to be 
positive when the flow is cyclonic, i.e., when the contributions 
of shear and rotation to the flow vorticity have the same sign. 
With our choice of axes, this implies that S = -2S xy , where 
S ij = l/2(diUi j + djiiij) is the usual deformation tensor. The 
system is closed either with the usual continuity equation sup- 
plemented by a polytropic equation of state, or, for simplicity, 
through an incompressibility assumption (V • w = 0). 

The relevant global time-scales of the problem are the shear 
time-scale t s = \S \, the viscous one t v = d 2 /v (d is the gap 
in the experiment, or the shearwise size of the shearing sheet 
box), and the rotation time-scale related to the Coriolis force 
tQ = (2Q)~' ; they relate to the advection term, the viscous term, 
and the Coriolis force term, respectively. Correlatively, the flow 
is described by two dimensionless numbers, the Reynolds num- 
ber 



Re = ty/t s = \S \d 2 /v, 
and the rotation number 

tfn = sgn(S)f J /f n = 2Q/S. 



(2) 



(3) 



For Keplerian flows, R& = -4/3. More generally, if one as- 
sumes that the large scale rotation of an astrophysical disk fol- 
lows a power-law, Q(r) oc r~ q , one locally has Rn = -2/q in 
the disk. 
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Note that our Reynolds number is defined on the outer 
scales, and not on the turbulent ones, such as, e.g., the Taylor 
microscale. Large values (~ 10 4 ) of this number are involved 
in the problem investigated here; the correlative numerical re- 
quirements are discussed in section l4~4l 

2.2. Energy budget 

As the global energy budget plays some role in the discussion 
of the results, it is rederived here. In the following equations, 
the bracket notation refers to a volume average of the brack- 
eted quantity. The averaging volume is the simulation one, and 
shearing-sheet boundary conditions are assumed in the deriva- 
tion, for definiteness. For the kinetic energy in the streamwise 
and shearwise directions, one finds: 



-M)=S (Rn + l){w x w y ) 



{ — — — ) + v(w A Aw A >, 
N p ox 1 



(4) 



^(y) =-SR n { w x w y ) 



'( —)+ V(WyA\Vy). 

v p oy 1 



(5) 



Instead of the vertical equation, it is more instructive to write 
down the total kinetic energy equation: 



where 



d iw 2 \ 

J t { — )=S(W X Wy)- 



e = vJ]<(Vw,) 2 > 



(6) 



(7) 



is the usual energy injection rate of turbulence cascade argu- 
ments 2 . In this last equation the incompressibility condition 
and the boundary conditions have been used in the reexpres- 
sion of the pressure term, and an integration by part has been 
performed on the viscous term (a constant kinematic viscosity 
v is assumed). 

In statistical steady-state, Eq. reduces to, 

S(\V X Wy)-£. (8) 

As pointed out bv lBalbus et alJ Jl996l) . the fact that e > 
implies that in steady state, the shear rate and the Reynolds 
stress responsible for radial transport have identical signs. 
This result has a direct physical interpretation: the imposed 
shear prevents the flow to be in global thermodynamic equilib- 
rium. Nevertheless, the flow tries to restore this global equilib- 
rium by radially transporting momentum through the turbulent 
Reynolds stress from regions of larger momentum to regions of 
lower momentum, consistently with Eq. (|8j. 

2 Because the rate of energy transfer in scale is constant in a 
Kolmogorov-like argument, the injection rate is directly related to the 
small-scale dissipation rate. 



Note finally that, in Eqs. @ and Q, the pressure-velocity 
correlation terms cannot be neglected, as they are of the order 
of the cascade energy injection term e. This is almost unavoid- 
able, as pressure is the only force that can provide for the ac- 
celeration of fluid particles in turbulent motions. As a matter 
of fact, the energy budget of any particular velocity compo- 
nent depends critically on the behavior of the velocity -pressure 
corre lations, which are notoriously difficult to model ( Speziale 
l!99lh . Ignoring this term in the analysis of the energetics there- 
fore leads to dubious or erroneous conclusions. 

2.3. Linear stability limits 

Surprisingly enough, the question of the linear stability lim- 
its of the simple rotating shear flows considered here is not 
completely solved to date. Focusing for the time being on 
purely streamwise-independent perturbations, i nstability with 
respect to local perturbations follows when llPedlevI Il969t 
lEehlanc & Camborll997tlsir^ & .TaccminkOOO) 



Rn(Rn + 1) < 0, 



(9) 



or, equivalently, -1 < Rn < 0. 

In plane Couette flows, it has been proven that R^ = 
is the correct cyclonic marginal stability limit for non 
streamwi se-invariant pe rturbations as well, at all Reynolds 
numbers ( Romanov Il973h . No such generic proof exists at the 
anticy clonic marginal stability limit (R n = -1). However, var- 
ious linear and nonlinea r numerical investigations suggest that 
this is indeed the ca se (Cam boii et alJll994tlKomminaho et alJ 
Il99rl iBech & Anderssonll997h . These results belong to plane 
Couette flows with rigid boundary conditions in the shearwise 
direction, but tend to prove that a local criterion captures the 
correct stabi l ity lim it, a s observed, e.g., in t he simulations of 
iBalbus et alJ i ll 9961) and lHawlev etaDJl999t) . 

The physics behind Eq. (|9ll can be captured by a disp laced 
particle argument ( Tritt on & Daviesl l 1981: Tritton 1992). This 
argument is reproduced in AppendixlAlfor the reader's conve- 
nience. Note that Eq. © is identical to Rayleigh's specific an- 
gular momentum criterion for the centrifugal instability, as the 
usual epicyclic frequency reads k 2 = S z Rq(1 + Rq). However, 
in the plane shear flow limit of cylindrical flows, the con- 
cept of specific angular momentum used in the derivation of 
Rayleigh's criterion no longer has meaning, so that one must 
follow a different route, as done here. Note also that, conse- 
quently, the Rayleigh criterion for the centrifugal instability in 
the inertial frame can also be understood from the action of the 
Coriolis force in the rotating frame (a somewhat surprising, al- 
though not new conclusion), as the displaced particle argument 
of Appendix|X|is readily extended to cylindrical flows. 

2.4. Subcritical transition in rotating plane Couette 
flows: a summary of relevant experimental results 

In the laboratory, non-rotating plane Couette flows undergo a 
subcritical transition to turbulence at Re 1500. The transi- 
tion Reynolds number steeply increases if a stabilizing rotation 
and/or a curvature is superimposed on the flow. The concep- 
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tually cleanest way to add rotation to a plane Couette flow is 
to place a plane Couette apparatus on a rotating table. Also, 
by considering a Taylor-Couette apparatus with varying gap 
width and independently rotating cylinders, one obtains a flow 
in which both rotation and curvature effects can be studied, and 
which reduces to a rotating plane Couette flow in the narrow 
gap limit. For a more complete discussion of the distinction and 
characterization of rotation and curvature in Taylor-Couette ex- 
periments, and_of_fliejeJate^_eM)erimental data, the reader is 
referred to Long aretti & Dauchotlll2005h . 

For the range of parameters studied to date in the exper- 
iments, it turns out that rotation and curvature effects on the 
transition Reynolds number are superposed in an mostly addi- 
tive way, so that both plane Couette flows and Taylor-Couette 
flows can in principle be used to characterize the effect of rota- 
tion. Concerning cyclonic flo ws, the only direc t ly rele vant data 
have been collected by Tillmark & Alfredsson ( 1996) with the 
help of a plane Couette flow apparatus placed on a rotating ta- 
ble. For anticyclonic flows, the only availa ble experiments are 
those of Richard and coworkers (lRichardl|200lt iRichard et alJ 
12001 h . who used a Taylor-Couette apparatus. The range of ro- 
tation number Rn explored in these experiments is — 0.1 for 

cyclonic rotation, and —1.6 1 for anticyclonic rotation. The 

data are shown on Fig. [2] 

The important point to note here is the steep dependence of 
the transition Reynolds number with the "distance" to marginal 
stability, with a typical slope \AR g \/\AR a \ ~ 10 4 — 10 5 . 




3. Numerical codes, strategy, and results 

In the present work, we are concerned with rotating, unstrat- 
ified uniform shear flows. Periodic boundary conditions hold 
in the direction of the flow (x axis) and the "vertical" direc- 
tion (z axis), and either shearing sheet or rigid boundary con- 
ditions are applied in the direction of the shear (y direction). 
The vertical axis is also the axis of rotation of the flow. The 
shearing shee t boun dary conditions are described in detail by 
lHawlev et al.l d 19951) . Shearing sheet flows thus modelled can 
be viewed as a local approximation of disk flows, while the use 
of mixed rigid-periodic boundary conditions is appropriate to 
numerically represent the rotating plane Couette flows of lab- 
oratory experiments, as routinely done in the fluid mechanics 
community. 

3.1. Numerical codes 

Two different 3D codes have been written for the present 
work: a finite difference compressible code, similar to ZEUS 
dStone & Normanl|l992l) . but restricted to the cartesian geom- 
etry, and rigid-periodic or shearing sheet boundary conditions; 
and a 3D incompressible Fourier code, in cartesian geometry, 
and implementing only the shearing sheet boundary conditions. 
An explicit kinematic viscosity term is added in both codes, 
upon which the Reynolds number is defined. Both codes were 
parallelized using the Message Passing Interface. 

The shearing sheet boundary conditions induce some 
changes with respect to a standard Fourier code. As a mat- 
ter of fact, while we were developing this code, the work 



R n - R n 

Fig. 2. Data on the Reynolds number of subcritical transition to tur- 
bulence as a function of the rotation number R n , measured from the 
appropriate marginal stability limit R „ (see text). Top panel: cyclonic 
plane Couette flow (data from Tillmark & Alfredsson 1996). Bottom 
panel: anticyclonic Taylor-Couette flow (data from iRichardEoOlh . 
The anticyclonic data are more difficult to collect, and consequently 
noisier. 



bv lUmurhan & Reeevl (Eo04) appeared, which implements the 
same technique. Therefore, our description of the required 
changes will be brief, and we refer the reader to this recent 
paper for details. 

To get effective periodic boundary conditions on the 3 axes, 
one needs to write Eq. Q in the sheared frame defined by: 



= t 



x + S -y-t 

y 
z. 



(10) 

(11) 

(12) 
(13) 



In this shearing frame, Eq. Q (supplemented by the incom- 
pressibility condition) becomes: 
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dw - 

— ¥ W ■ VlO 

dt' 



V ■ w 



V6n 

P 

+vAw 
0. 



■ 2Q w x e y + (2Q + S )w y e x 



(14) 
(15) 



in which V = d^e^ + (dy - S t'd x >)e y ' + d z >e z - and A = V ■ V. 

Since the shearing box is a periodic box in the shearing 
frame, this last formulation of the Navier-Stokes equation can 
be written in 3D-Fourier Space. Defining 



one finally obtains: 



dw .. 

¥ iu ■ w ® w 

dt' 



fi — k - S tk x e y 



8n 

-ifi 2Q w x e v + (2Q + S)w x e a 

P 

-v/i 2 w 



(16) 



fi ■ w = 



(17) 
(18) 



These are the equations actually used in our spectral code. 
The nonlinear term is computed using the 2/3 dealiasing rule 
with a pseudo-spectral method (see e.g Pevret 2002 for a de- 
scription of this point) and each time-step is evaluated using 
a 4 th order Runge Kutta Scheme. One should note that a k- 
wave in the sheared frame actually appears as a yu(f)-wave in 
the steady frame. Then, as time goes on in the simulation, the k- 
grid describes higher spatial frequency in the steady frame and 
consequently, the large scales are not computed anymore. Since 
nonlinear coupling limits the shearing of any wave-number, 
a remap procedure is periodically applied all along the sim- 
ulati on, and pr events to loose information on the large scale 3 
(Rogallo 198 l j). This kind of algorithm has been extensively 
described bv llJmu rhan & Regev (2004) using a 2D spectral 
code and the reader should refer to this publication for tech- 
nical details on the remap procedure. 

The choice of these two codes was made first for purposes 
of comparison with previous work, and secondly to allow us 
to cross-check the potential limitations of one code against the 
other; e.g., the shearing sheet boundary conditions and sheared 
spatial basis Eq. (I16> have their own limitations, as the sheared 
basis forms a complete basis for shearing sheet boundary con- 
ditions, but only for these conditions. 

The three codes were tested in a variety of ways. The first 
test was to reprodu ce the non-rotating pla ne Couette flow be- 
havior computed bv lHamilton et al.M 1995b . This was done both 
with our finite difference code, and with David Clarke's version 
of ZEUS3D, for comparison purposes. We checked the non- 
linear transition mechanism was well reproduced, with the cor- 
responding Reynolds number and aspect ratio, and that the two 
codes gave completely consistent results. Then, the shearing 
sheet boundary conditions were tested using these two finite 
difference codes and the Fourier code. We have verified that 
mean turbulent quantities (e.g., mean energy, mean transport, 
velocity maxima and minima) and critical Reynolds number 

3 We thank Achim Wirth for pointing out this reference to us. 



were statistically the same using the different codes, for dif- 
ferent rotation numbers, either cyclonic or anticyclonic. This 
consistency holds over the 10 5 - 10 6 time steps of our simula- 
tions. 

3.2. Initial conditions and numerical strategy 

The experimental results recalled in sec tion l2~4l suggest that a 
steep dependence of the transition Reynolds number with the 
rotation number may be the cause of the difficulty to find such 
a transition in the previously published shearing sheet numer- 
ical simulations. Accordingly, one of the major aims of this 
investigation is to quantify the effect of the simulation resolu- 
tion on the determination of the transition Reynolds number as 
a function of R&. 

Now, one of the characteristic features of the sub- 
critical transition to turbulence is an observed spread in 
transition Reynolds numbers, depending on the choice of 
initial conditions, and a correlative large spread in tur- 
bulence life-t imes. This has been docu mented both ex- 
perimentallv foarbvshire & Mullinl 1 19951) and numerically 
( Faisst & Eckhardt 2004) in pipe flows, and guides to some 
extent our choice of initial conditions and our numerical proce- 
dure. Indeed, turbulent life-times typically vary from fast decay 
(survival for less than one hundred dynamical times) to long or 
indefinite survival (several thousands of dynamical times, with 
a clear divergence at finite Reynolds number) over several or- 
ders of magnitude of variation of the initial condition ampli- 
tude, _butforjess_dian50% of variations of Reynolds number 
(see Faisst & Eckhardt 2004, Figs. 2 and 7). 

It is reasonable to assume that this qualitative behavior is 
generic. Consequently, we have chosen once and for all, fixed, 
high amplitude initial conditions, to make our numerical runs 
more directly comparable to one another upon variations of 
Reynolds numbers. Furthermore, we consider that turbulence 
is long-lived if it is not observed to decay for 100 or 200 
shear times (depending on the runs). This choice is a com- 
promise between computational time constraints, and accuracy 
in the determination of the transition Reynolds number of in- 
definitely self-sustained turbulence. In practice, simulations are 
performed in a cubic box (the impact of this choice is discussed 
in the next section, to some extent). The flow is adimension- 
alized with the only dimensional quantities introduced in the 
problem: S and d, where d is the simulation box size (or equiv- 
alently, by choosing \S\ — 1 and d — 1). The initial conditions 
used for all our simulation are a random 3D excitation of the 10 
largest Fourier modes, with rms fluctuations in velocity of or- 
der unity in our chosen units. Other shapes of initial conditions 
were tested such as white noise (all scales excited randomly) 
or introducing large scale vortices in various directions with a 
small superimposed noise. This produces no significant differ- 
ence once the flow is relaxed (f > 20 t s ). 

The numerical strategy adopted is then rather straightfor- 
ward: choosing a code, a resolution, a boundary condition (for 
the finite difference code) and a Reynolds number, at fixed ini- 
tial conditions, the flow evolution is computed starting from the 
marginal stability limit in rotation number Rq and evolving the 
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rotation number by (small) fixed steps every 100 or 200 shear 
times. According to the preceding discussion, this allows us to 
reduce at maximum the number of runs and the run time needed 
to observe systematic trends in the numerical results. 

In this section, only shearing sheet boundary conditions are 
used. We have also checked that the time required to dissi- 
pate the turbulent energy of the flow assuming energy injec- 
tion is stopped (deduced from the e term in Eq. (6)) is smaller 
than 100f. s ; this constraint is always satisfied by a large mar- 
gin in all our runs, implying that the deviations from laminar 
motion that we observe are self-sustained (i.e., we do not ob- 
serve them because their dissipation time exceeds the run time). 
Actually, once turbulence is lost in our simulations, the energy 
in the velocity fluctuations always decreases rather fast, as can 
be checked on Fig.[3]for cyclonic flows. The same property is 
found for anticyclonic flows, see section l3~4l 

We conclude this section on our choice of the Mach num- 
ber (Ma = L V S /c s ) for our simulations with the compressible 
Zeus-like code. The type of motions we are considering in these 
simulations reach at most a small fraction of the boundaries rel- 
ative velocity (normalized to unity in this work). We found that 
a sound speed also normalized to unity was a good compromise 
between limiting the effects of compressibility (which eventu- 
ally makes the turbulence compressible and largely different 
in character when the Mach number is too large), and the im- 
pact of the sound speed on the CFL condition. Also, this value 
mimics the real role of compressibility in a vertically stratified 
accretion disk. Consequently, Ma = 1 is imposed in all our 
compressible simulations. 

3.3. Numerical results: cyclonic flows 

On the cyclonic side, simulations are performed while main- 
taining the rotation number Rq constant during 100 shear 
times; then the rotation number is increased by steps of 0.02, 
starting from the marginal cyclonic point Rq = 0. An exam- 
ple global output of such a simulation is plotted on Fig.[3]for 
Re= 12000. The relaminarization point is easily found since the 
transition between the turbulent to laminar state is quite abrupt 
(at t = 1150 on Fig. 0. We define the last turbulent point 
as the last rotation rate for which turbulence is sustained for 
100 f s . For our example simulation, we find that the last tur- 
bulent point at Re= 12000 and 64 3 resolution with our Fourier 
code is Rq = 0.2. 

Using this kind of simulation, we plot the last turbulent 
points in the (Re, Rq) space, for different resolutions and/or 
codes on Fig. |4] Turbulence is found on the cyclonic side at 
least up to Rq = 0.3, i.e significantly away from the marginal 
stability point. 

Note that turbulence is maintained with certainty (with our 
adopted criteria) at any given point, but, due to the sampling 
made in the explored Reynolds number, turbulence may also 
be maintained at a somewhat lower Reynolds number (i.e just 
below the last turbulent point on Fig. [3}. This can be true down 
to the previously tested Reynolds number, for which turbulence 
is not maintained at the considered rotation rate. In conclu- 
sion, the real transition Reynolds Rg curve in the (R&, Re) plane 
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Fig. 3. Example of the time evolution of a 64 3 (Re= 12000) cyclonic 
flow run as computed by our Fourier code. The turbulent energy, trans- 
port and dissipation rate are the quantities involved in Eq. {6}. The 
dissipation time follows from the turbulent energy and the dissipation 
rate. The bottom panel displays the evolution of the rotation number 
that is imposed in the course of the simulation. 



should be found somewhat below (but not far from) the last 
turbulent point curve determined here. This remark is more im- 
portant for anticyclonic flows, for which precise quantitative 
results are needed. 

Except for a systematic shift between the results obtained 
with the Fourier code and the ZEUS-like one, the results seem 
to be independent of the resolution. The numerically minded 
reader may ask how one can reach such high Reynolds numbers 
with such relatively small resolutions. This point is addressed 
in section l4~4l 

An important issue is to quantify transport in subcritical tur- 
bulent flows. The phenomenological arguments of Lon garettil 
( 2002) suggest that (v x v y ) oc 1 /Rg in subcritical flows, and that 
the turbulent transport in a given flow with specified (Re, Rq) 
numbers depends only on Rq through Rg (see section 14. 1> 4 . 
Consequently, we have used all our simulations at a given R& 
to obtain the least noisy evaluation of (v x v y ). Then, with the 
help of Fig. [3] one finds a transition Reynolds number Rg for 



4 The same result follows if one assumes that in the fully turbulent 
state, the torque oc Re 2 , as predicted in Kolmogor ov turb ulence, and 
observed in experiments (see, e.g.. lDubrulle et al. 2005a). The argu- 
ment of section lPl allows us to recover this result from more generic 
physical principles. 
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Fig. 4. Transition Reynolds number Rg as a function of the rotation 
number R n , with different resolutions and codes for shearing-sheet 
boundary conditions (cyclonic rotation). All points were obtained us- 
ing our Fourier code except those labelled FD (finite difference) which 
use our ZEUS-like code. 
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Fig. 5. Mean transport as a function of the transition Reynolds num- 
ber for cyclonic rotation (normalized by S 2 d 2 ). 



any given Rn, which allows us to plot the mean turbulent trans- 
port (v x v y ) as a function of the transition Reynolds number in 
Fig-0 This was done only from the data of our Fourier code for 
self-consistency, but using both the 32 3 and 64 3 resolution runs, 
as they produced the same results, and as the use of a larger data 
set improves the statistics. The resulting relation reads 



(v x Vy) 



5.5 



Rg - 1250 



(Sd) 2 



(19) 



The presence of an additive constant in the denominator of this 
expression is a clear indication of the influence of the linear in- 
stability close to the marginal stability limit; indeed, transport 
in the supercritical region is significantly enhanced with respect 
to the subcritical region (see, e.g., Fig. 16 in iDubrulle et alJ 
2005a, and explanations therein). For large critical Reynolds 
number (i.e., far enough from the marginal stability boundary, 
e.g., Rg > 15000), (v x v v ) 5.5/Rg is a good approximation. 



3.4. Numerical results: anticyclonic flows 

The strategy adopted in simulations of anticyclonic flows is 
similar to the cyclonic side. Starting at Rn = -1.0, the rotation 
number is increased in steps of -0.004 and each step lasts 200 
shear times to allow for flow relaxation. A typical run is shown 
on Fig.|6] computed with our 3D Fourier Code at Re = 12000. 
One should note that the flow fluctuations have higher ampli- 
tudes on the anticyclonic side than on the cyclonic side; this is 
why we have reduced the rotation number steps and increased 
the relaxation time in anticyclonic runs. Consistently, The last 
turbulent point is defined here as the last rotation rate for which 
turbulence is sustained for 200f s . On the example Fig. |6j we 
find the last turbulent point for Re = 12000 at R n = - 1 .024. 

As for cyclonic rotation, the last turbulent points for anti- 
cyclonic rotation are plotted on Fig. and the mean transport 
on Fig. [F] Error bars are added on Fig. [H]to help assessing the 
significance of the various fits performed, as they will be used 
later on. On the lower bound of these bars, turbulence is not 
maintained with certainty whereas the contrary is true for at 
least 200 shear times at the upper bound. Therefore, the actual 
transition Reynolds number is bracketed by the error bar. 

Recalling that Rn = -2/q with D.(r) oc r~ q and that 
Rn — -1 corresponds to a constant specific angular momen- 
tum distribution in cylindrical flows, the largest rotation num- 
ber reached here (Rn = -1.032) corresponds to q — 1.94; 
this is quite consis tent with the results shown on Fig. 1 of 
lHawlev et alJ Jl999l) . except for the crucial fact that the reso- 
lution and Reynolds number dependence are now quantified. 
The reason why such high Reynolds numbers are accessible 
with our relatively low resolutions is discussed in section l4~4l 
For the time being, let us comment a bit further on the informa- 
tion encoded in Fig. [7] which shows that Reynolds number and 
resolution are different, albeit related control parameters. We 
will focus on the Fourier code data for definiteness. Consider 
the 32 3 data, for example. For \Rn\ < 1.016, the transition 
Reynolds number agrees with the one found at higher resolu- 
tion. However, increasing the Reynolds number above ~ 6000 
produces a loss of turbulence at the same rotation number inde- 
pendently of the Reynolds number, whereas this is not true at 
higher resolutions. This implies that the physics is not faithfully 
represented at this resolution for Re > 6000 and Rn > 1 .016. 

This is the most important point to note here: two different 
regimes of transition from turbulent to laminar are displayed 
in this figure. The first (corresponding to the various fitting 
curves) is the correct, resolution independent and Reynolds de- 
pendent transition. The second (apparent as the various ver- 
tically aligned points at a given resolution) is an incorrect, 
Reynolds independent and resolution limited transition. Note 
that the points belonging to both this vertical line and the 
laminar-turbulent line are still resolved, though, as shown in 
section 14.4.21 The meaning of the behavior displayed in Fig.0 
is further discussed in section |4~TI and its implications in sec- 
tions E2| and E3I 

Comparing Figs. 0] and we remark that the dependence 
of the transition Reynolds number Rg on the "distance" to 
marginal stability in rotation number \Rn -R%\is considerably 
stiffer on the anticyclonic side than on the cyclonic one. This 
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Fig. 7. Transition Reynolds number Rg as a function of the Rotation 
number Rq, and related analytical fits, with different resolutions and 
codes for shearing-sheet boundary conditions (anticyclonic side). All 
plots were computed using our Fourier code except FD (finite differ- 
ence) which uses our ZEUS-like code. Note that the x-axis is inverted 
with respect to Fig. [2] Symbols along the fitted lines correspond to 
resolved simulations; vertically aligned symbols indicate the limiting 
rotation number that can be reached at a given resolution and mostly 
correspond to unresolved simulations. For the sake of clarity, symbols 
which sit on top of each other have been slightly displaced along the 
Rq axis; this is indicated by the arrows and the related values of Rq. 



Fig. 6. Time evolution of a 64 3 Re= 12000 anticyclonic flow as com- 
puted by our Fourier code. Panel description is identical to Fig.|3| 



has important implications that will be discussed in the next 
section. Conversely, the turbulent momentum transport is very 
similar to the one found for the cyclonic side 5 , as shown on 
Fig-Gfl 



5.5 



-JSd) 2 



(20) 



The constant in the denominator differs from the one found 
on the cyclonic side. This reflects the difference of transition 
Reynolds number at the two marginal stability limits. For large 
enough Reynolds number, one find (v x v y ) - 5.5 /Rg, which 
corresponds to the asymptotic relation found on the cyclonic 
side (see section P~TI for a discussion of the possible origin of 
this behavior). This indicates that this relation is very robust 
for subcritical flows, far enough from the supercritical transi- 
tion limit. 

4. Discussion 

Our results are at variance wi th b oth the point of view ad- 
vocated bv lBalbusetai1lll996l> and lHawlev et alJ i 19991) (ab- 
senc e of subcritic a l turb ulence), and Richard & Zahn (1999) 
and iHersant et al.l d2005l) (efficient transport due to subcriti- 
cal turbulence). This is further investigated in this section. We 

5 Fig. [5] is noisier than its cyclonic counterpart. This is a conse- 
quence of the larger turbulent fluctuations observed in anticyclonic 
flows. Longer integrations time-scale would have been required to im- 
prove the statistics. 
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Fig. 8. Mean transport as a function of critical reynolds number on 
the anticyclonic side (normalized by S 2 d 2 ). 



shall first present some phenomenological background mate- 
rial which helps to understand the physical origin and mean- 
ing of the results presented in the previous section. Then, we 
shall respectively discuss the implications of our results for 
Keplerian flows (section l4~2l . the stabilizing role of the Coriolis 
force in subcritical flows (section l4~3l . and the relation between 
Reynolds number and resolution (section l4~4l : these last two 
items have been highly controversial in the past decade. Section 
I4.4I also discusses the relation of the se results with the scale- 
invariance argument of Balbus (2004). Finally the influence of 
the nature of the adopted boundary conditions and aspect ratio 
on our results is the object of section l4~5l as well as and their 
relation to fluid dynamics experiments. Note also that the dis- 
cussion of the boundary conditions helps quantifying possible 
biases introduced by the sheering sheet boundary conditions 
with respect to actual disk physics. The reader interested only 
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in the astrophysical implications of our results may focus on 
sectionlllll 



4.1. Some aspects of subcritical turbulence 
phenomenology: 

The phenomenology of subcritical turbul ence has been dis - 
cussed in Longaretti (2002) and Lonaaretti & Dauchot (2005). 
Some directly relevant aspects for our present purpose are pre- 
sented here (and clarified where needed). 

Turbulent transport is often quantified in terms of a turbu- 
lent viscosity. As this description has been criticized in the past, 
a brief discussion of its use here might be useful. First, note 
that, in scale-free systems such as the ones studied here (the 
only scale present being the simulation box size), one can al- 
ways assume that 



(v x v y ) = v,S, 



(21) 



as this only amounts to defining a turbulent viscosity v, such 
that this relation is satisfied. In any case, as the source of tur- 
bulence is the shear, the Reynolds stress (v^Vy) must be some 
function of the shear S , which cancels when the shear cancels. 

Now, v, has the dimension of a length times a velocity, so 
that one must therefore have, in our simulations, 



v, = aSd , 



(22) 



as S d and d are the only dimensional quantities with the right 
dimensionality introduced in the problem. 

ff is a Shakura-Sunyaev-like parameter. It is a dimension- 
less quantity, and can therefore only depend on the dimension- 
less quantities 6 characterizing the problem at hand, namely the 
Reynolds number Re and the rotation number /fo (i.e., the shear 
dependence of a can only appear through the ratios of the shear 
time scale to the viscous and the rotation time scales): 



a = a(Re,Rn). 



(23) 



The results of section [3] suggest that, quite remarkably, a 
depends only on Rq through the transition Reynolds Rg, and 
not (or little) on Re, in subcritical flows. The ori gin of this 
behav ior can be understood in the following way JLong aretti 
120021) . 

A sheared flow is out of global thermodynamic al equilib- 
rium, and tries to restore this equilibrium by transporting mo- 
mentum across the shear. A subcritical flow has only two means 
at its disposal to achieve this purpose: laminar and turbulent 
transport. It will tend to choose the most efficient one under 
any given set of conditions 7 , i.e. at given Re and Rq. The sub- 
critical turbulent transport will exceed the laminar one when 



6 Actually, in principle, a depends also on the aspect ratio of the 
simulation, and on the nature of the boundary conditions. As these are 
not varied in the results discussed on the basis of the phenomenology 
described here, this dependency is ignored for simplicity. 

7 Note that this does not imply that the momentum transport is ab- 
solutely maximized. 



v, > v. Right at the laminar-turbulent threshold, Re 
v, ~ v. This implies that 

v 1 
a - . 

Sd 2 Rg 



Rg and 



(24) 



Now, what does happen at Reynolds numbers Re larger 
than the transition Reynolds number Rg ? To answer this ques- 
tion, it is useful to have in mind some idealized, qualitative 
picture of the situation in wave-number space. Such a picture 
is proposed in Fig. [9] and constitutes a reasonable working hy- 
pothesis. It is reasonably well-supported by our current knowl- 
edge of the plane Couette flow turbulent self-sustaining pro- 
cess and of inertial spectra, as well as by the spectral analysis 
of some of our simulations presented and discussed in section 
I4A2I 
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Fig. 9. Proposed sketch of the idealized energy spectrum in a turbu- 
lent shear flow. Arrows indicate the energy flow through mode cou- 
pling. The length-scales l M and lj correspond to the top of the inertial 
range (assumed identical to the bottom of the self-sustaining range for 
simplicity) and the top of the dissipation range. Scales are assumed to 
be normalized to the box simulation size d, and anisotropy is ignored 
in this sketch (see text for details). 



In this picture, the large scales are occupied by the self- 
sustaining mechanism. All scales in this domain are expected to 
be coherent in phase, and interactions between large and small 
scales occur both ways 8 . The intermediate range is the inertial 
range of turbulence; scales have no phase coherence, energy 
cascades to smaller scales at a constant rate, provided by the 
self-sustaining mechanism (as part of the mode coupling tak- 
ing place in the self-sustaining mechanism range of scales oc- 
curs with the inertial range). The smallest range represents the 
viscous dissipation scales. The existence of the self-sustaining 
process scales, their properties, and their influence on the iner- 
tial range (energy input and anisotropy) is the distinctive fea- 



8 This is the case in particular for the non-rotating plane Couette 
self-sustaining mechanism 1 Waleffe 1997). 
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ture of shear turbulence with respect to the more commonly 
known and studied forced isotropic turbulence. 

In such a picture, increasing the Reynolds number almost 
exclusively results in an increase of the inertial range, which is 
essentially vanishing at the transition Reynolds number. This 
should have little effect on the turbulent transport (whereas, on 
the contrary, the laminar transport becomes smaller and smaller 
when increasing the Reynolds number). 

Indeed, we have first checked that this is case in non- 
rotating Couette plane flows, where the self-sustaining mecha- 
nism is identified (Hamilton et al. 1995): the transport is almost 
completely determined dominated by the mechanism rolls and 
streaks. Furthermore, in our simulations, we have computed the 
contribution of each length scale to the total transport (v v v v ). 
First one should note that in Fourier space (in ID for simplic- 
ity): 

N-l 

(v x v y ) = J] v x (k n )v;(k n ) (25) 
«=o 

Therefore, the contribution to mean transport of the wavelength 
k„ is found to be 25l(v v (£„)v*(A:„)), since k„ and k^^ n represents 
the same physical wavenumber. This simple result can be used 
in 3D by averaging the transport over 2 directions (in physical 
space) and by computing the Fourier transform in the remain- 
ing direction ; this procedure is sufficient for our purpose here. 
The resulting cumulative Fourier sum, starting from n — is 
illustrated on an example on Fig. ^| to quantify which scales 
dominate the transport. In this example, the resulting spectral 
analysis is plotted in the x direction for a 128 3 anticyclonic 
flow with Rq = -1.024 and Re = 20000, showing that more 
than 99% of the transport comes from scales larger than 1/10 
of the box size ; this range corresponds to the length-scales of 
the self-sustaining process (see section l4.4.2i . Similar results 
are found for spectral analyzes in the y and z directions, consis- 
tently with the picture discussed here. This is expected anyway 
if the inertial spectrum is kolmogorovian, as confirmed from 
the spectral analysis of section 9 14.4.21 

There are two loose ends in this discussion. First, hysteresis 
is usually experimentally observed in subcritical transitions to 
turbulence: the measured transition Reynolds number is higher 
when moving "up" from the laminar to turbulent states than 
when moving "down" from the turbulent to laminar ones. This 
suggests that the laminar-turbulent boundary is separated by 
some sort of barrier in the appropriate phase-space (defined, 
e.g., by the amplitudes and phases of the Fourier modes). This 
(along with the fact that the arguments developed here apply 
only in order of magnitude) may well explain the existence of 
the constant of order 5 that one finds in Eqs. dl9> and ( 1201 with 
respect to Eq. ( I24l >. Secondly, the arguments presented here 
ignore the existence of marginal stability thresholds. This, as 
pointed out in sections l3~3l and l3~4l may explain the presence 
of the constant at the denominator of these relations, as the 

9 The true nature of the inertial spectrum might be affected by the 
anisotropy generated by the the shear and the Coriolis force, but these 
anisotropies must become negligible at small scale, due to the shorter 
and shorter eddy turnover time. 
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Fig. 10. Example of the cumulative sum contribution of each scale 
length to the mean transport, starting from k x = in the x direction: 
99% of the transport comes from k x < 10 scales. From a 128 3 simula- 
tion, Re=20000, R n = -1.020. 



equivalent global subcritical transition Reynolds number that 
one can define in the supercritical regime is orders of magni- 
tude smaller than in the subcritical regime. 

To conclude, let us point out the relation of this picture with 
the numerical results presented in Fig. [7] The fact that higher 
resolutions are required to faithfully represent the physics at 
higher rotation numbers indicates that the ratio cL/Im increases 
with rotation number. Indeed, if the resolution is too low, so 
that the relative scale l M /d is not resolved, the energy transfer 
loop represented on Fig. [5] cannot take place, and turbulence 
is not self-sustained. Furthermore, at the transition Reynolds 
number, the inertial spectrum is nearly inexistent, as pointed 
out above, and Im ~ Id- Consequently, the most critical scale 
ratio in this problem is expected not to be the Kolmogorov one, 
but the self-sustaining mechanism one (if //m)- 

4.2. Implications for Keplerian flows 

Actual disks are vertically stratified, whereas stratification is ig- 
nored in our experiments. Stratification provides us with a local 
macroscopic scale (the disk scale height H). With appropriate 
provisos related to the possible stabilizing or destabilizing role 
of stratification 10 , one can tentatively identify this scale height 
with our simulation box size: H — d. This assumption is made 
throughout this section. In the same way, the Shakura-Sunyaev 
ass parameter is defined such that v, = a S sc s H a ss QH 2 . 
Eq. J22i then implies that a S s = 2a7l^nl a (the last equal- 



10 If stratification is destabilizing, the momentum transport induced 
by the resulting convective motions is in the wrong direction, as re- 
called in the introduction, and must be counterbalanced by another 
process; ignoring stratification in this case therefore makes life eas- 
ier for this other process (here, subcritical turbulence). If stratification 
is stabilizing, this also most likely results in an increased difficulty 
in finding the transition to turbulence, and a related increase in the 
transition Reynolds number. These arguments suggest that ignoring 
the dynamical stratification altogether maximizes the overall outwards 
transport in our problem. 
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Table 1. Extrapolated transition Reynolds numbers, values of 
a, and required simulations resolution, for keplerian flows, un- 
der various assumptions (see text for details) 





exponential 


power-law 


cyclonic 


linear 


Rg 


1.3 x 10 26 


1.1 x 10 8 


2 x 10 7 


1.8 x 10 6 


a 


n/a 


5 x 1(T 10 


2.6 x 1(T 7 


3.1 x 10- 6 


(d/6) 3 


n/a 


7000 3 


3000 3 


900 3 



ity holds within a factor of order unity for the rotation number 
range of interest in this work). 

Using the numerical results shown on Figs. [7] and [8] one 
can deduce a few properties of Keplerian flow subcritical shear 
turbulence, based on various conservative extrapolations of our 
numerical data. First, the transition Reynolds number Rg de- 
pendence on the rotation number R& is well-fitted by a power 
or an exponential law. Using these laws, one can get a first set 
of estimates of the transition Reynolds number for keplerian- 
like flows (R n = -4/3): Rg = 1.1 x 10 10 and Rg = 1.3 x 10 26 , 
respectively. The last estimate leads to the absence of subcrit- 
ical turbulence in accretion disks whereas the first one allows 
for its existence 1 1 . Secondly, let us note that, for both cyclonic- 
ity, the Coriolis force induces a steeper and steeper increase of 
the transition Reynolds number when moving away from the 
marginal stability boundary. This suggests that one can find a 
lower bound for Rg by linearly extrapolating the power law fit 
beyond the last known point (Rn = -1.032). One find this way 
Rg mm = 1.8 x 10 6 . As a final hypothesis, one may envision 
that the Rg(Rn) relation would be more or less symmetric with 
respect to Rq = if there were no supercritical domain. This 
would explain why the actual relation of Fig. is so steep: in 
this picture, the system tries to reach back as fast as possible 
the high values of transition Reynolds number expected from 
this hypothetical symmetry, after which the Reynolds depen- 
dence with rotation number would be much less steep. Under 
this assumption the expected transition Reynolds number for 
keplerian flows would be Rg = 2. x 10 7 (a power-law fit of the 
cyclonic data has been used in this extrapolation). 

This information is summarized in tabled along with the 
corresponding values of a, obtained from the asymptotic rela- 
tion a = (v x v y ) = 5.5 /Rg found for cyclonic and anticyclonic 
flows in the previous section. The last line shows the resolution 
required to successfully simulate keplerian flow turbulence, for 
the various Reynolds numbers (see subsection l4.4.1> . One sees 
that even the most optimistic a bound (a max = 3.1 x 10~ 6 ), 
obtained with the linear extrapolation, is substantially smaller 
than the values req uired in astrophysical accretion disks (as 
summarized, e.g., in lPapaloizou & LirJl995l) . Note finally that, 
even without any extrapolation, our results exclude subcritical 
turbulent transport at the a 3.10~ 4 level. 



1 1 We assume that accretion disk Reynolds numbers lie between 10 10 
and 10 15 for definiteness. The Reynolds number definition used in this 
evaluation is Re = SH 2 /v where H is the local disk scale height, con- 
sistently with the H = d identification made earlier. 



4.3. Role of the Coriolis force in uniform shear flows 

Two different but related issues have been raised in the liter- 
ature concerning the role of the Coriolis force in subcritical 
systems. 

[ for linearly stable flows, Bal bus et alJ l l 19961) point 
out that the Coriolis force plays a conflicting role in Eqs. @ 
and (|5j- More precisely, they make the following point: as 
S(v x v y ) > for turbulence to exist (see section l2~2l . the terms 
in which the shear S has been factored out in these equa- 
tions have opposite signs for linearly stable flows, while they 
have the same sign for linearly unstable flows (note that this 
is true independently of the flow cyclonicity). They conclude 
from this that a stabilizing rotation prevents turbulence to show 
up in subcritical shear flows, except possibly in the vicinity 
of marginal stability. So mewhat relatedly, t he recent Reynolds 
stress-closure model o f lOgilvid J2003I) and lGaraud & Qgilvid 
(2005) predicts relaminarization for large enough deviations 
from the marginal stability limit. In particular, for the au- 
thors' standard choice of parameters, it predicts relaminariza- 
tion for Rq ~ 0.2 for cyclonic rotation. H owever, as can be 
seen on Fig. |4| both thelBalbus et alJ i ll 9961) argument and the 
iGaraud & Ogilviel J2005I) result conflict with our simulations: 
subcritical turbulence is maintained away from marginal stabil- 
ity on the cyclonic side, at least up to Rq 0.3. Note that we 
could have pushed the search for transition to turbulence be- 
yond what is shown on this graph, especially by using higher 
resolutions, but did not do it due to computer resources limita- 
tions. As discussed in the next subsection, t he absence of turbu- 
lenc e in the keplerian flo w simulations of Balb us et alJ (Ti 996) 
and lHawlev et alJ Jl999) is a problem of resolution. 

The second issue relates to the asymmetry between cy- 
clonic and anticyclonic rotation. The stress-closure model just 
mentioned depends on the rotation number only through the 
combination Rn{Rn + 1) which implies a symmetry with re- 
spect to Rq = -1/2. This symmetry is clearly violated by our 
numerical results (compare Figs. @] and Q, a point which re- 
quires some comments. 

First, note that the linearized Navier-Stokes equation Q ex- 
hibits this symmetry for perturbations with vanishing pressure 
variation (Sn = 0). In this case, the linearized equation can be 
written: 

~dt = s ' y !h + s ' ( (Rn + )WyCx ~ RnWxe y) + yAw (26) 

The cyclonic-anticyclonic symmetry appears when exchanging 
the x and y directions. Indeed, upon the following change of 
variables: 





= -Rn-l, 






< 


= Wy, 


e', 


= Cy. 




= w x , 




= e* 


w' 


= w z , 


e', 


= e,, 



so that 

w' = w' x e' x + w' y e'y + w' z e' z 
= w, 
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the form of Eq. i26\ should be invariant, which is indeed the 

case: 

?T = s ' y % + s ' ( (R ' n + 1)w > ,e '- r " + vAw '' (27) 

This symmetry can also be extended to compressible motions 
by adding 6n'{x,y, z) = 6n(y,x,z) to the list of change of vari- 
ables. 

Because the perturbations defining the linear stability 
limit also exhibit this symmetry (Appendix [XJ, it has often 
been assumed in closure-stress models in the past. However, 
this is not a symmetry o f the full Navie r -Stokes equation 
JSneziale & Mhuirisl 1 19891 ISneziald fl99ll ISalhi & Gambol 
199?|, nor of the V • w — equation). This is also appar- 
ent in a direct inspection of the structure of simulated turbu- 
lent flows. The Rn = 0, wall-bounded turbulent flows con- 
tain la rge streamwise rolls living for about a hundred shear 
times (Hamilt on et alJlT995l) . We have also found rolls more 
or less aligned in the streamwise direction in our Rq = 
shearing sheet simulations, although we did not try to precisely 
quantify their survival time. Furthermore, at the anticyclonic 
marginal stability limit (Rn = -1), we did observe sheared 
shearwise rolls (i.e rolls in y direction) in our simulations, as 
one might expect from the symmetry of the linearized Navier- 
Stokes equation. The anticyclonic roll survival time is observed 
to be rather short compared to their cyclonic counterpart, as 
they are tilted by the shear and loose their coherence in a few 
shear times at most. This roll lifetime is the main difference 
we found between the cyclonic and anticyclonic side. This is 
related to the fact that a streamwise roll does not reduce the 
shear on the anticyclonic subcritical domain (in opposition to 
the cyclonic one). 

In any case, we have found turbulence away from the 
marginal stability limit in cyclonic flows, and the symmetry 
with respect to R& = -1/2 is violated both in our simulations, 
and in supercritically r otating shear flow turbulence (see, e.g., 
ISalhi & Cambon 1997 and references the rein). Th i s mak e the 
predictions of the stress- closure model of Oeilvie (2003) and 
Idaraud & Ogilviel |2005 ) unreliable in both subcritical and su- 
percritical flows. 

4.4. Resolution, effective Reynolds number and scale 
invariance 

The results of section l3~3*l and l3~3l invol ve fairly high Reynolds 
numbers, and one might ask if our simulations are resolved 
enough in these regimes. This question has a priori two dif- 
ferent aspects, as one can guess from Fig. [9] resolving the self- 
sustaining process smallest relative scale cI/Im, and resolving 
the relative dissipation scale d/lj. 

For the problem considered in this paper, resolving the first 
scale is a sine qua non condition: if it is not satisfied, turbulence 
does not show up, independently of the simulation Reynolds 
number, because the required scale coupling shown on Fig. |9| 
for the self-sustaining process to exist cannot take place. This 
shows up in Fig.0as the vertical transition limit from turbulent 



to laminar that we obtained for any given resolution, for large 
enough Reynolds numbers. 

Resolving the dissipation scale is important to ascertain that 
direct numerical simulations such as the ones performed here 
are not biased by (the presence or absence of) numerical dis- 
sipation, and this issue is often raised in the fluid mechanics 
literature. For the time being, we note that, at the transition 
Reynolds number, the inertial domain should be non-existent 
or extremely reduced, so that Im — h an d both resolution re- 
quirements should be directly related (this point, used in sec- 
tion is justified in section 14.4.21 . We can therefore con- 
sider that the "effective Reynolds number" Re^ of our simu- 
lations is the largest transition Reynolds number Rg correctly 
determined at a given resolution 12 , as discussed in section IX?! 

Note that this effective Reynolds number is problem- 
dependent: the self-sustaining process qualitative and quanti- 
tative characteristics both depend on the considered problem; 
furthermore, in simulations of isotropic turbulence, the self- 
sustaining process is absent, and replaced by a forced ampli- 
tude of the largest Fourier modes, so that the effective Reynolds 
number in this case is the one related to the dissipation scale. 

Let us now examine the two requirements mentioned above 
in more detail. 

4.4.1. Resolving the self-sustaining process 

First, we would like to qualitatively comment on the difference 
of resolution requirements between cyclonic and anticyclonic 
flows. 

As discussed in section I4.5I the nature of the shearwise 
boundary condition has apparently only a small influence 
on the results; this is exemplified by the similar transition 
Reynolds numbers found in our simulations and in experi- 
ments on rotating shear flows (see Fig. 114-b . This suggests that 
at least some of the characteristics of the self-sustaining pro- 
cess of non-rotating plane Couette flows are relevant here. At 
the cyclonic marginal stability limit, th is self-sustaining pro- 
cess has a tim e-scale fssp ~ 1005 _I JHamilton et alJ ll995: 
IWaleffell9 97). The requirement that, at the transition Reynolds 
number, the viscous time scale at scale Im exceeds fssp reads 
l 2 M /v > 1005"', i.e., l„/d < (100 /Rs) 1/2 ~ 1/4 for Rg ~ 
1500 dLongaretfi & Dauchotl 120051) . This probably explains 
why the resolution requirement is so low on the cyclonic side. 
Conversely, we have mentioned at the end of the previous sub- 
section that rolls (which are an apparently ubiquitous ingredi- 
ent in subcritical turbulence) do not survive more than a few 
shear times in anticyclonic flows. Therefore, the anticyclonic 
self-sustaining process time-scale cannot exceed a few shear 
times as well, whatever its nature. The same reasoning as the 
one exposed above leads to Im/<J < a few (1/Rg) 1 ^ 2 ~ a few 
xl/70, an already much more demanding constraint. It is ob- 
viously related to the larger transition Reynolds number found 
at the anticyclonic marginal stability, compared to the cyclonic 
one. 

12 With all the provisos discussed in section l3~2l about the role of 
the choice of the initial conditions and turbulence minimal survival 
life-time. 
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(d/6) 3 




r 


32 3 


6000 


1.71 


64 3 


12000 


1.48 


128 3 


38000 


1.46 



Table 2. Resolution, effective Reynolds number and y factor for the 
Fourier code on the anticyclonic side 



As mentioned several times already, the self-sustaining pro- 
cess is identified and understood only at the cyclonic marginal 
stability limit in wall-bounded Couette flows. Consequently, it 
is difficult to explain why the resolution demand grows so much 
faster with rotation number "distance" to marginal stability for 
anticyclonic flows than for cyclonic ones. However, we spec- 
ulate that this is connected to the fact the rotation time scale 
is only a fraction of S~ l for cyclonic flows, whereas it always 
exceeds S~ x for anticyclonic ones. 

Next, let us try to quantify the resolution that would be 
needed to successfully simulate keplerian flows. The phe- 
nome nology of subcritical turbulence developed by Lon garettil 
J2002I) predicts that d/l M ~ Rg 1 ' 2 and (v x v y ) oc 1 /Rg. This phe- 
nomenology implicitly assumes that the relevant time-scale of 
the self-sustaining process is ~ S , so that it would need to 
be modified to be applied to cyclonic flows, but it should be 
adequate for anticyclonic ones, with appropriate modifications. 
In particular, we have already pointed out in section 14.11 that 
the last relation needs to be amended into (v x v y ) oc l/(Rg - R c ) 
with R c si 3000 on the anticyclonic side. This suggests that 



6 
d 



(Rg-R c ) 1 ' 2 



(28) 



is the appropriately generalized scale relation (6 being the 
smallest scale accessible to the simulation, i.e., the resolution). 
Table |2] gives the values of y and Re^ for the three different 
resolutions of our anticyclonic simulations. 

Although the statistics is a little poor to draw firm con- 
clusions, it appears that y is nearly constant compared to the 
variations in both resolution and transition Reynolds number, 
and our simulations are therefore consistent with Eq. (1281 . The 
resolution needed to simulate keplerian flows has been com- 
puted based on the estimate Eq. d28i . with y — 1 .5 (the R c cor- 
rection has little influence on these estimates). The results are 
shown in tabled For comparison purposes, note that the largest 
turbulence simulation ever performed was 4000 3 , but was not 
run for hundreds or thousands of dynamical times. Although 
the results gathered here are probably only indicative, as they 
are based on guess work, they strongly suggest that simulating 
subcritical turbulence in keplerian flows is beyond present day 
computer capabilities, and support the idea that the subcritical 
keplerian flows simulations performe d to date were l imited by 
numerical resolution, as suggested by Longaretti ( 2002). 

4.4.2. Resolving the dissipation scale 

In statistically steady turbulence, the dissipation scale can be 
defined from the balance between input and dissipation de- 
scribed by Eq. <|HJ. The energy input is provided by S{v x v y ). 



The Fourier analysis of this quantity is shown on Fig. EH an d 
is dominated by the large scales. Conversely, the Fourier con- 
tent of e, Eq. 0, is dominated by the small scales (large k), 
comparable to the dissipation scale, as illustrated below. 

Resolving the dissipation scale is important with Fourier 
codes in order to prevent energy accumulation at the smallest 
scales, which may bias the results, or lead to code crash 13 . 

The general definition of the dissipation wavelength kj fol- 
lows from the evaluation of Eq. Q in Fourier space: 



Jrkd 




k 2 E(k)dk 



(29) 



where it is assumed that E(k) is cut-off at kj, either abruptly, 
or through so me modelling of the dissipation range (see e.g. 
lLesieuJl990l) . 

In simulations of homogeneous and isotropic turbulence, 
the energy input is imposed from the outside: the amplitude of 
the largest Fourier mode is held fixed, and Fig. |9] reduces to 
the inertial and dissipation range. In this context, the inertial 
spectrum reduces to the Kolmogorov spectrum given by: 



E K (k) = C K e 2/3 k- 5/ \ 



(30) 



where the Kolmogorov constant Ck - 1. Cutting off this spec- 
trum at wavelength kd and injecting it in the definition Eq. J29I 
leads to the well-known expression of the Kolmogorov wave 
number, k/c = (e/v 3 ) 1 ^ 4 . The related Kolmogorov scale (inverse 
of the wave number) is a largely used estimate of the dissipa- 
tion scale. 

In the fluid mechanics community one often requires that 
the Kolmogorov scale be resolved, even if the considered turbu- 
lence is not isotropic and ho mogeneous, as, e.g., in shear flow 
turbulence (see, e.g., |PumiH rT996). However, in our simula- 
tions, the observed spectrum is substantially different from the 
Kolmogorov one, especially at the transition Reynolds num- 
ber (see top panel of Fig. [121. Indeed, at the turbulent-laminar 
transition, one does not expect nor observe the presence of an 
inertial domain in the spectrum. One may therefore ask what 
relation the Kolmogorov scale bears to the dissipation scale of 
the problem. 

Consider, e.g., the 32 3 and 64 3 energy spectra obtained at 
a Reynolds number Re = 6000 and a rotation number set to 
-1.016. These spectra are shown on the top panel of Fig. 1121 
The concordance of the spectra at both 32 3 and 64 3 resolu- 
tions indicates that the dissipation scale in the 32 3 simula- 
tion is resolved (this is consistent with the shape of the spec- 
trum at the smallest 32 3 resolved scales, much steeper than 
Kolmogorov). It appears that the largest distance to marginal 
stability \Rq + 1 1 reliably accessible at a given resolution on the 
laminar-turbulent transition (as checked by higher resolution 
simulations) corresponds to the various vertical line of transi- 
tion displayed on Fig.[7]for this resolution. In other words, the 
Re = 6000, R n = -1.016 point at 32 2 , and the Re = 12000, 
Rq - -1.024 at 64 3 , are resolved. This feature makes us con- 
fident that the transition point determined at 128 3 is the correct 



13 One may also include an hyper- viscosity term to prevent code 
crash, but this turned out not to be necessary. 
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N 


R = k d /k K 


Rg 


32 


1.23 


6000 


64 


1.73 


12000 


128 


2.66 


35000 



Table 3. Resolution, dissipation to Kolmogorov wave number ratio, 
and corresponding transition Reynolds number (see text for details). 



one, although we did not cross-check it at 256 3 , due to the lim- 
itations in the available computational resources. 

We have thus determined the largest transition Reynolds 
number where the dissipation scale is confidently resolved in 
these anticyclonic runs at the various resolutions we have used 
(32 3 , 64 3 and 128 3 ). In other words, we know the effective 
dissipation scale of these simulations, as it must be compara- 
ble to the largest wave number available in the simulation 14 : 
kj 3 l ^ 2 nN/d, where N is the resolution. Furthermore, we 
can compute the Kolmogorov wave number k% for these runs, 
asv = Re/Sd 2 , and as e follows from Eq. (|8j and the transport 
(e.g., with the help of the transport/transition-Reynolds-number 
correlation displayed on Fig.|8j. The resulting ratio R = kj/kic 
is given in table[3] 

Although the values of the ratio R quoted in table are 
of order unity, a systematic trend seems to appear, indicating 
that resolving the Kolmogorov wave number is possibly not 
the relevant concept at the transition Reynolds number, as it is 
not stringent enough; nevertheless, the required resolution de- 
rived from the Kolmogorov wave number is apparently semi- 
quantitatively correct, at least for the rotation numbers explored 
here. Of course when going to Reynolds numbers well in ex- 
cess of Rg at a given Rq, the Kolmogorov wave number should 
always give the right estimate of the dissipation scale, as the 
inertial range becomes more and more prominent in the overall 
spectrum. 

To conclude this aspect of the discussion, we note that 
both the non-kolmogorovian shape of the spectrum at transi- 
tion and the relatively small values of e at the various transition 
Reynolds numbers used here, most probably combine in the 
end to explain why we can reach rather large Reynolds num- 
bers at rather moderate resolutions. 

In order to have a better grasp on which scales contribute 
most to the dissipation, we have computed a quantity, Td(k), 
defined by 

Xk pk pk 
dkj dkyl dk z {k\ + k 2 + k 2 )E(k x ,k v ,k z ). (31) 
k J-k J-k 

Comparing with Eq. i29\ . it appears that Td(k) represents 
the fraction of dissipation due to scales \k x \ < k, \k y \ < k and 
\k z \ < k. This quantity is plotted on Fig. ^2 with the 64 3 simu- 
lation spectrum. It appears that more than 95% of the total dis- 
sipation is due to k < l/2k m . dx (i.e., the 32 3 resolution). Also, 
comparing Fig. II II with the top panel of Fig.ll2lindicates that 

14 This expression corrects a misprint in Pumii 1 1996) for the diag- 
onal of a cube in Fourier space; although this largest wave number is 
resolved only in discrete directions, this definition is adopted here for 
ease of comparison with this earlier work. 
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Fig. 11. Cumulative mean dissipation spectrum for a 64 3 simulation 
at Re=6000 for R a = -1.016. 

most of the dissipation comes from the part of the spectrum 
which is steeper than the Kolmogorov spectrum, as one would 
expect. 

It is also instructive to examine the spectral behavior at 
Reynolds number larger than the transition Reynolds number, 
as shown on Figll2l 

This figure displays energy spectra of the velocity devi- 
ation from the laminar flow. The rotation number is fixed at 
Rn = -1.016 for all spectra, and they have been averaged over 
a 200 shear time period to reduce the noise. From top to bottom, 
the Reynolds number is 6000, 12000 and 20000 respectively. 
The 32 3 simulations are resolved only in the top panel, while 
the 64 3 simulations should be resolved in the top two panels. 
Comparing the second panel with the first reveals a couple of 
interesting points: 

- The 64 3 simulation shows an extension of the spectrum, 
compatible with a small inertial range (this is difficult to 
ascertain because of the remaining noise in the simulation), 
while still resolving at least the top of the dissipation range, 
but marginally so. 

- The 32 3 simulation begins to significantly deviate from the 
64 3 simulation, although the trend is similar. 

The third panel also displays a fairly relevant piece of in- 
formation. The 64 3 simulation shows both the self-sustaining 
mechanism scales and the inertial spectra. However, the dissi- 
pation scale does not seem to be resolved. This is not unex- 
pected, since increasing the Reynolds number necessarily in- 
creases the inertial spectrum, and therefore decreases the dissi- 
pation scale. Apparently, the dissipation scale is probably not 
far off the resolved scales, so that the simulation nevertheless 
does not noticeably deviate from the expected behavior. But 
note that the 32 3 simulation is clearly strongly unresolved, with 
energy accumulating in the small scales in order for a statisti- 
cally steady equilibrium to be achieved: indeed, as this simula- 
tion resolves the self-sustaining mechanism scale, turbulence is 
present; however, as the smallest resolved scale is significantly 
larger than the dissipation scale, the spectrum must be strongly 
deformed to achieve a dissipation which is consistent with the 
energy input due to the turbulence self-sustaining mechanism. 
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Fig. 12. Energy spectra (of the velocity deviation from the laminar 
flow), for two different resolutions (32 2 and 64 3 ). The rotation number 
is R n = -1.016 in all cases. Top panel: Re = 6000. Middle panel: Re = 
12000. Bottom panel: Re = 20000. At this resolution only the top 
panel simulations are resolved, as expected. See text for discussion. 



These simulations illustrate that if the dissipation scale is 
not resolved, the simulated flow does not necessarily relami- 
narize, but the deformation of both the amplitude and shape of 
the spectrum most likely results in, e.g., unreliable estimates 
of the turbulent transport. In particular, the reliability of finite 
difference simulations where no viscous term is explicitly in- 
cluded in the code is unclear 15 . 

On the other hand, this point is related to the fact that 
the numerical dissipation in a Fourier code is extremely weak, 
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We did not further investigate this question here. 
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Fig, 13. Energy budget for a 64 3 Re = 20000 run with our Fourier 
code. Each plot represent a term in Eq. JSJ. The numerical dissipation 
is normalized by the total dissipation d(w 2 /2)/dt - S{w x w y }. We find 
that numerical dissipation is about 1% of the total dissipation. 



so that the deformation of the spectrum may be quite large. 
To compute the numerical dissipation explicitly, we have esti- 
mated its effect on the turbulent energy budget. 

We plot an example of such an energy budget on Fig. \^\ 
where all the terms in Eq. are evaluated, so that the remain- 
ing difference measures the code dissipation. One should note 
that these plots are integrated over 2 shear times, so that they 
include the numerical dissipation due to the desaliazing pro- 
cedure (done at each time loop) and losses from the remap- 
ping procedure (done each shear time). The presented result is 
generic: for all our simulations, numerical dissipation is found 
to be at most a few percent of the total dissipation. 

In summary, we have tried as much as possible to ensure 
that our determination of the transition Reynolds number and 
turbulent transport do not suffer from lack of resolution of 
the dissipation scale. Note also that the results of the Fourier 
and finite difference codes are consistent with each other. This 
makes us confident that our simulations faithfully represent the 
relevant physics, down to and including the dissipation scale, 
within the relevant (Re, Rn) domain determined at each resolu- 
tion on Fig. [7] 

4.4.3. Shearing sheet simulations and scale 
invariance 

Recently, Balbus| d2004 has argued that the scale invariance 
of the inviscid form of the Navier-Stokes equation used here 
makes any small scale solution exist at large scales as well. 
This argument seems to imply that simulations of the kind per- 
formed here shoul d not be resolution limited. However , neithe r 
the simulations of Bal bus et all Jl99fj|) . iHawlev et all |l999), 
the ones performed here, nor a real disk, are scale invariant. In 
shearing sheet simulations, the box size defines a scale; in a real 
disk, the disk scale height does. Furthermore, we point out that 
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the mechanism analyzed bv lWaleffel Jl997l) . whose qualitative 
and semi-quantitative relevance to the present work has been 
pointed out hereabove, is somewhat insensitive to the nature of 
the imposed boundary condition. Along with the results found 
in this paper, this suggests that only a scale rather than a spe- 
cific boundary condition needs to be imposed for statistically 
stationary turbulence to show up in numerical simulations, as 
exemplified in section [3~4l Finally, the role of an increasingly 
dominant Coriolis force is not to define another scale, which it 
cannot, but to modify the relative range of scales that are re- 
quired for turbulence to exist (most likely because of its more 
and more stringent time-scale requirement), so that numerical 
resolution does play an important role in subcritical turbulence 
detection, as can be seen from Fig.0 
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Fig. 14. Rg(Ro) plot from experimental data iTillmark & Alfred sson 
1996, crosses), and our numerical simulations using 64 3 Fourier code 
(circles) and 64 3 finite difference code (triangles) with cubic box and 
shearing sheet boundary conditions. 



4.5. Boundary conditions and aspect ratio 

Assessing the role of boundary conditions on the existence and 
properties of subcritical turbulence is an important question, 
since real accretion disk boundary conditions are not repro- 
ducible in experimental flows. However, the resolution demand 
in the local shearing box is already so large for a keplerian flow 
that a global simulation of a subcritical keplerian disk flow is 
totally out of reach. The best we can do is to compare numer- 
ical experiments with shearing sheet and rigid/periodic bound- 
ary conditions with one another, and with experimental results. 
This is the object of this section. 

Before doing so, let us point out some important differences 
between the two types of boundary conditions: 

- In the semi-lagrangian variables defined by Eqs. dl0> . the 
only difference is that the velocity deviation from the lami- 
nar flow cancels on the rigid boundary in the rigid/periodic 
case, while it is periodic in the shear direction (as in the oth- 
ers) in the shearing sheet case. This results in a suppression 
of the boundary layer in the shearing sheet case. 

- Characteristic sizes are the same in both cases. However, 
while for rigid/periodic conditions, structures are forced to 
remain more or less stationary with respect to the walls on 
average, this is not the case with shearing sheet boundary 
conditions, where structures can move at random through 
the boundary. As a consequence, a long-lasting mean flow 
distortion is apparent with rigid/periodic boundary condi- 
tions (due to the matching of turbulently enhanced trans- 
port with the viscous one in the boundary layer), while 
in shearing-sheet simulations, although such a distortion is 
usually locally found at any given time, it averages out over 
time, due to its random localization. 

- This relates to a profound difference between accretion 
disks and actual experiments. In the latter, the flow pro- 
file adjusts to the imposed boundary condition through a 
pressure redistribution, and a stationary state is reached. In 
the former, this cannot take place, and the disk is never sta- 
tionary, due to the resulting turbulent transport of mass and 
angular momentum. 

In spite of these differences , we shall nevertheless argue 
that the choice of boundary conditions has only a limited im- 
pact on some of our qualitative and semi-quantitative results. 



This suggests that the underlying mechanisms are reasonably 
closely related in both settings, although much more work than 
what has been possible to do here is required to ascertain this 
conclusion. 

4.5.1 . Cyclonic rotation 

Fig. 1141 displa ys a comparison o f our numerical results with 
the lTillmark & Alfredssonl fl996) data, in the range of rotation 
number where these data were collected. 

The agreement between the two is fair, with the Fourier 
code results being sensibly more compatible with the data than 
the finite difference code ones, at the larger rotation numbers. 
This follows because, at the same "resolution", a Fourier code 
is physically more resolved than a Finite difference code. Note 
also that some 128 3 simulations were performed using the 
Fourier code and the same transition thresholds were found 
as for the 64 3 simulations. This supports the idea that the 64 3 
Fourier code results are not resolution limited. 

We have also made a few runs using rigid (shearwise direc- 
tion) and periodic (other directions) boundary conditions with 
our ZEUS-like code. At each rotation number, we made a few 
tries with different Reynolds numbers to locate the transition 
threshold. Each run was computed from the same initial con- 
dition for 400 shear times with 80 x 80 x 40 grid points and 
a L x — 1 .757T Ly = 1 L z = 1 .2n aspect ratio box (correspond- 
ing to the "minimal flow unit" aspect ratio, i.e. the smallest 
box in which turbulenc e can be sustain ed with these bound- 
ary conditions: see lHamilton et alll 19951 for details). The error 
bars upper bounds correspond to the lower Reynolds for which 
turbulence is found and the lower bound the higher Reynolds 
number for which turbulence is lost. The numerical data are 
shown on Fig. [21 the error bars reflect our poor sampling, not 
intrinsic fluctuations in the transition Reynolds number. These 
data are fitted by a linear law: 

Rg = 1400 + 4 x l0 5 R n , (32) 

the slope of which is 15 times steeper than the one found from 
the experimental data. 

This dramatic difference in transition Reynolds number 
with respect to the experimental and shearing sheet results is 
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Fig. 15. Re c as a function of R n for cyclonic rotating plane couette 
flow. 



in fact controlled by the choice of the simulation box aspect 
ratio. For example, let us choose a longer box in the z direction 
(i.e L x = \.15nL y = 1 L z — 2 An). With such a choice, turbu- 
lence is sustained at Rn = 0.01 and Re = 2400, much closer 
to the expected transition R eynolds of Fig.fTHthan what is pre- 
dicted by Fig. ^] Finally, iKomminaho et alJ J1996I) . using a 
very elongated simulation box in the flow direction, found tran- 
sition right at the experimentaly determined Reynolds number 
(Rg = 3000, R n = 0.06). 

These result show the important role of aspect ratio in 
subcritical turbulence simulations with rigid/periodic bound- 
ary conditions. Apparently, the use of shearing sheet boundary 
conditions relaxes this constraint. This is reasonable since the 
shearing sheet box allows more freedom than rigid boundary 
conditions. In actual experiments, the aspect ratio is not an is- 
sue since usually very large L x /L y and L z /l y are used, so that 
the turbulence coherence length can freely adjust itself in these 
directions. 

These results also indirectly suggest that the turbulence sat- 
uration mechanism is not strongly affected by the use of shear- 
ing sheet boundary conditions. One would nevertheless expect 
that the reduction of the shear in the middle of the flow, due 
to the mean velocity profile modification which occurs with 
rigid/periodic boundary conditions, produces a reduced turbu- 
lent transport. This is indeed the case: e.g., the turbulent trans- 
port at marginal stability (R n = 0) is (v c v v ) 2 x lO^(Sd) 2 
for the rigid/periodic boundary conditions 16 , while one has 
(v x v y ) - QA(Sd) 2 throughout the flow with the shearing sheet 
boundary conditions, although the transition Reynolds number 
is the same in both instances. These features most probably 
find a natural explanation if the turbulence amplitude satura- 
tion mechanism is primarily controlled by the system nonlin- 
earities, and not by the mean profile deformation. 

4.5.2. Anticyclonic rotation 

The comparison of our numerical results with Richa^feOJHl) 
data is shown on Fig. 1161 

The discrepancy between the experimental and numerical 
data is striking, especially at the light of the remarkable con- 

16 This is measured in the middle of the flow where the turbulent 
transport is maximized, and viscous transport negligible. 
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Fig. 16. Rg(Rn) plot from experimental data on Taylor-Couette flows 
l Rich arcfefJOlL crosses), and the various numerical simulations results 
and related fits shown on Fig.Q 



sistency observed for cyclonically rotating flows. In particular, 
the increase in transition Reynolds is considerably steeper in 
the numerical data than in the experimental ones. Note however 
that the numerical and experimental data seem to give the same 
transition Reynolds number at the marginal stability boundary. 

lLongaretti & Daucholl 12005) have argued that the flow cur- 
vature plays no role in the anticyclonic flow data of iRicharcl 
(2001), so that the origin of the large discrepancy between the 
numerical and experimental results must be found elsewhere 17 . 
In this respect, note that experimental secondary flow distor- 
tions are much more likely to induce a linear instability some- 
where in the flow on the anticyclonic side as on the cyclonic 
one. Indeed, recall that the stability limit is defined by Eq. 
Consider the cyclonic marginal stability limit first (Rn = 0), 
and assume that one moves away from it by imposing a small 
change in rotation SO.. The required change in shear profile SS 
to locally achieve 26Sl/(S(y) + SS) < is large: SS ~ S(y) is 
needed. Conversely, at the anticyclonic marginal stability limit 
(Rn = -1, i.e., S = 2D.), upon a small change SO. of the ro- 
tation rate, a change SS - 2SQ. <sc S suffices to locally make 
2Q./S > -1 and produce a linear instability somewhere in the 
system. This argument shows that the presence of secondary 
flows, such as Ekmann's circulation, can easily make the flow 
more unstable than it would be in its absence in anticyclonic 
flows, whereas this is much more difficult to achieve in cy- 
clonic ones. This may easily explain the discrepancy between 
numerical and experimental results shown on Fig. [HO while the 
agreement is remarkable at the marginal stability boundary. 

5. Summary and conclusions 

The central results of this paper are displayed on Figs.@][5]0 
|S]and|5] and their significance and implications are discussed 
in sections EHl l3~4l |4~T1 IO R31 andg3] The main implica- 
tions of these results are summarized in the abstract. In the 
course of the discussion, we have found that a number incor- 



17 In any case, the flow curvature always increases the transition 
Reynolds number, so that including curvature in the analysis of this 
problem can only make it worse, not better. 
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rect statements have been made in the literature, most notably 
concerning the existence and importance of subcritical turbu- 
lence in presence of a dynamically significant Coriolis force. 
We have also found that resolution is a key issue for subcritical 
anticyclonically rotating flows (including keplerian ones), and 
have quantified the relation between resolution, rotation and 
Reynolds number. In relation to this, we believe that the ques- 
tion of resolution of the dissipation scale is not emphasized 
enough in the astrophysics literature, and the potential effects 
of this problem are most probably underestimated. 

Our simulations do not faithfully represent a real disk: nei- 
ther vertical stratification, nor, more critically, realistic vertical 
boundary conditions have been implemented. A real (hydrody- 
namic disk) moves either in vacuum, or, more probably, in a hot 
corona. In both cases, one expects the real vertical boundary 
condition in the disk to be (nearly) stress-free. We have made 
some very preliminary simulations of ah stratified disk-corona 
system to test this idea, where most of the inertia lies in the 
disk. Although a strong numerical mixing of the corona and 
the disk at the interface prevents us to evolve the system for 
a long time (50t s max), no significant difference in the over- 
all dynamics of the disk did show up. However this problem 
probably requires more careful investigations to validate this 
conclusion. 

Overall, the outcome of this investigation still leaves us 
with the issue of transport unresolved in MHD-inactive flows 
(and possibly in some MHD-active ones), and we will briefly 
comment the various ways out of this conundrum. 

We first note that an efficient enough local instability should 
lead to a large enough turbulent transport, because the tran- 
sition to fully developed turbulence usually occurs to signifi- 
cantly lower Reynolds numbers in these systems than the ones 
found here. This is true, e.g., in rotating shear flows of the type 
considered here, in the linearly unstable regime. However, no 
such instability has yet been found in hydrodynamical keple- 
rian flows, either stratified or not, as discussed in the introduc- 
tion. It remains to be seen whether another such instability can 
operate in hydrodynamic disks, but the list of potential driving 
agents has by now significantly been narrowed. 

In what concerns the YSO disks dead-zone in particular, 
it may be that the disk stirri ng due to the M RI above and be- 
low the dead-zone itself Fleming & Stand 12003) might pro- 
vide enough transport in the end if it excite s large eno ugh large 
scale 2D disturbances of the right type ( Ioan nou & Kak ouris 
l200ll) in the disk. However, this option remains to be worked 
out in detail. 

It has often been noted that transport in disks may not be 
due t o turbulence but to waves (see, e.g., Papaloi zou & Linl 
1 1 995l for an introduction to the subject). R ecent results on the 
existe nce of vortices in stratified disks ( Barr anco & Marcusl 
120051) and on the coupling of wav es to vortices resu lting in 
efficient transport in 2D dynamics jRodo et alJl2005l and ref- 
erences therein) support this idea. 



Appendix A: Displaced particle analysis for 
rotating flows: 

The following line of argume nt closely follows 

iTritton & Daviesl Jl98lh and iTrittonl Jl992l) . Let us con- 
sider a rotating shear flow, whose dynamics is controlled by 
Eq. Q. As in section |2~T1 x is the direction of the flow, y the 
direction of the shear, and z the direction perpendicular to the 
x,y plane, in which the rotation ii is applied. The laminar 
equilibrium velocity Ul = (U(y), 0,0) generates a Coriolis 
force in the y direction of magnitude -2pQ.U (in algebraic 
value), which is balanced by the equilibrium generalized 
pressure gradient —dn/dy. 
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Fig. A.l. Sketch of the effect of the Coriolis force on the displaced 
particle (see text). 



Let us further consider two fluid "rods" of infinite extent 
in the streamwise direction x, and located at positions y\ and 
yi — y\ + Sy. The streamwise velocities of these rods are U\ 
and U2, respectively. Let us assume that one displaces the rod 
atji to location;^, without disturbing the pressure distribution. 
Although the total work of the Coriolis force vanishes, there is 
a net partial work due to the force component in the x direction 
which originates in the velocity v of this displacement in the 
y direction. Because of this partial work, the rod experiences 
a change of x momentum, and therefore of x velocity, which 
reads 

U[-U 1 = j 2Qvdt = 2Q6y, (A.l) 

so that the velocity U[ of the rod when it reaches location y^ 
differs from the equilibrium velocity U2, and correlatively, the 
x component of the Coriolis force acting on this displaced rod, 
-2pOt/J (in algebraic value) differs from the equilibrium one, 
-2pQXJ 7 (seeFig.lATV 

Consequently, the net result between the equilibrium pres- 
sure gradient and the Coriolis force will tend to restore the dis- 
placed rod to its equilibrium position 18 if U[ > U%, or displace 
it further if U[ < U^- From Eqs. JA. II and (0, one obtains 



U:-U 2 = 2Q.5y - —Sy = S(R n + l)Sy. (A.2) 
dy 

where S = -dU/dy is the shear. From this result, the net force 
(Coriolis and pressure) on the displaced rod reads 



2pQ.(U 2 -U[) = - P S 2 R n (Rn + l)6y, 



(A3) 



18 Consider the direction and magnitude of the Coriolis force and 
pressure force along the y axis to derive the effect of the inequality. 
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This shows that equilibrium is always restored when Rq > or 
Rn < -1 and destroyed otherwise (equality holds at marginal 
stability). This is the result quoted in section l2~Tl This result can 
also be directly derived from the linearized eulerian equation 
of motion with the use of spatially uniform perturbations of the 
pressure and the velocity. 
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